home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / GELS.FOR < prev    next >
Text File  |  1985-11-29  |  5KB  |  183 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE GELS
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH
  8. C           SYMMETRIC COEFFICIENT MATRIX UPPER TRIANGULAR PART OF WHICH
  9. C           IS ASSUMED TO BE STORED COLUMNWISE.
  10. C
  11. C        USAGE
  12. C           CALL GELS(R,A,M,N,EPS,IER,AUX)
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           R      - M BY N RIGHT HAND SIDE MATRIX.  (DESTROYED)
  16. C                    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
  17. C           A      - UPPER TRIANGULAR PART OF THE SYMMETRIC
  18. C                    M BY M COEFFICIENT MATRIX.  (DESTROYED)
  19. C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
  20. C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
  21. C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
  22. C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
  23. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  24. C                    IER=0  - NO ERROR,
  25. C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
  26. C                             PIVOT ELEMENT AT ANY ELIMINATION STEP
  27. C                             EQUAL TO 0,
  28. C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  29. C                             CANCE INDICATED AT ELIMINATION STEP K+1,
  30. C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
  31. C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
  32. C                             ABSOLUTELY GREATEST MAIN DIAGONAL
  33. C                             ELEMENT OF MATRIX A.
  34. C           AUX    - AN AUXILIARY STORAGE ARRAY WITH DIMENSION M-1.
  35. C
  36. C        REMARKS
  37. C           UPPER TRIANGULAR PART OF MATRIX A IS ASSUMED TO BE STORED
  38. C           COLUMNWISE IN M*(M+1)/2 SUCCESSIVE STORAGE LOCATIONS, RIGHT
  39. C           HAND SIDE MATRIX R COLUMNWISE IN N*M SUCCESSIVE STORAGE
  40. C           LOCATIONS. ON RETURN SOLUTION MATRIX R IS STORED COLUMNWISE
  41. C           TOO.
  42. C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
  43. C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
  44. C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
  45. C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
  46. C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
  47. C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
  48. C           GIVEN IN CASE M=1.
  49. C           ERROR PARAMETER IER=-1 DOES NOT NECESSARILY MEAN THAT
  50. C           MATRIX A IS SINGULAR, AS ONLY MAIN DIAGONAL ELEMENTS
  51. C           ARE USED AS PIVOT ELEMENTS. POSSIBLY SUBROUTINE GELG (WHICH
  52. C           WORKS WITH TOTAL PIVOTING) WOULD BE ABLE TO FIND A SOLUTION.
  53. C
  54. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  55. C           NONE
  56. C
  57. C        METHOD
  58. C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
  59. C           PIVOTING IN MAIN DIAGONAL, IN ORDER TO PRESERVE
  60. C           SYMMETRY IN REMAINING COEFFICIENT MATRICES.
  61. C
  62. C     ..................................................................
  63. C
  64.       SUBROUTINE GELS(R,A,M,N,EPS,IER,AUX)
  65. C
  66. C
  67.       DIMENSION A(1),R(1),AUX(1)
  68.       IF(M)24,24,1
  69. C
  70. C     SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT
  71.     1 IER=0
  72.       PIV=0.
  73.       L=0
  74.       DO 3 K=1,M
  75.       L=L+K
  76.       TB=ABS(A(L))
  77.       IF(TB-PIV)3,3,2
  78.     2 PIV=TB
  79.       I=L
  80.       J=K
  81.     3 CONTINUE
  82.       TOL=EPS*PIV
  83. C     MAIN DIAGONAL ELEMENT A(I)=A(J,J) IS FIRST PIVOT ELEMENT.
  84. C     PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
  85. C
  86. C
  87. C     START ELIMINATION LOOP
  88.       LST=0
  89.       NM=N*M
  90.       LEND=M-1
  91.       DO 18 K=1,M
  92. C
  93. C     TEST ON USEFULNESS OF SYMMETRIC ALGORITHM
  94.       IF(PIV)24,24,4
  95.     4 IF(IER)7,5,7
  96.     5 IF(PIV-TOL)6,6,7
  97.     6 IER=K-1
  98.     7 LT=J-K
  99.       LST=LST+K
  100. C
  101. C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
  102.       PIVI=1./A(I)
  103.       DO 8 L=K,NM,M
  104.       LL=L+LT
  105.       TB=PIVI*R(LL)
  106.       R(LL)=R(L)
  107.     8 R(L)=TB
  108. C
  109. C     IS ELIMINATION TERMINATED
  110.       IF(K-M)9,19,19
  111. C
  112. C     ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.
  113. C     ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX.
  114.     9 LR=LST+(LT*(K+J-1))/2
  115.       LL=LR
  116.       L=LST
  117.       DO 14 II=K,LEND
  118.       L=L+II
  119.       LL=LL+1
  120.       IF(L-LR)12,10,11
  121.    10 A(LL)=A(LST)
  122.       TB=A(L)
  123.       GO TO 13
  124.    11 LL=L+LT
  125.    12 TB=A(LL)
  126.       A(LL)=A(L)
  127.    13 AUX(II)=TB
  128.    14 A(L)=PIVI*TB
  129. C
  130. C     SAVE COLUMN INTERCHANGE INFORMATION
  131.       A(LST)=LT
  132. C
  133. C     ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT
  134.       PIV=0.
  135.       LLST=LST
  136.       LT=0
  137.       DO 18 II=K,LEND
  138.       PIVI=-AUX(II)
  139.       LL=LLST
  140.       LT=LT+1
  141.       DO 15 LLD=II,LEND
  142.       LL=LL+LLD
  143.       L=LL+LT
  144.    15 A(L)=A(L)+PIVI*A(LL)
  145.       LLST=LLST+II
  146.       LR=LLST+LT
  147.       TB=ABS(A(LR))
  148.       IF(TB-PIV)17,17,16
  149.    16 PIV=TB
  150.       I=LR
  151.       J=II+1
  152.    17 DO 18 LR=K,NM,M
  153.       LL=LR+LT
  154.    18 R(LL)=R(LL)+PIVI*R(LR)
  155. C     END OF ELIMINATION LOOP
  156. C
  157. C
  158. C     BACK SUBSTITUTION AND BACK INTERCHANGE
  159.    19 IF(LEND)24,23,20
  160.    20 II=M
  161.       DO 22 I=2,M
  162.       LST=LST-II
  163.       II=II-1
  164.       L=A(LST)+.5
  165.       DO 22 J=II,NM,M
  166.       TB=R(J)
  167.       LL=J
  168.       K=LST
  169.       DO 21 LT=II,LEND
  170.       LL=LL+1
  171.       K=K+LT
  172.    21 TB=TB-A(K)*R(LL)
  173.       K=J+L
  174.       R(J)=R(K)
  175.    22 R(K)=TB
  176.    23 RETURN
  177. C
  178. C
  179. C     ERROR RETURN
  180.    24 IER=-1
  181.       RETURN
  182.       END
  183.